Reverberating activity in a neural network with distributed signal transmission delays 
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It is known that an identical delay in all transmission lines can destabilize macroscopic stationarity 
of a neural network, causing oscillation or chaos. We analyze the collective dynamics of a network 
whose intra-transmission delays are distributed in time. Here, a neuron is modeled as a discrete- 
time threshold element that responds in an all-or-nothing manner to a linear sum of signals that 
arrive after delays assigned to individual transmission lines. Even though transmission delays are 
distributed in time, a whole network exhibits a single collective oscillation with a period close to 
the average transmission delay. The collective oscillation can not only be a simple alternation 
of the consecutive firing and resting, but also nontrivially sequenced series of firing and resting, 
reverberating in a certain period of time. Moreover, the system dynamics can be made quasiperiodic 
or chaotic by changing the distribution of delays. 
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I. INTRODUCTION 

A number of model neural networks exhibit collective 
oscillation. Their mechanisms of oscillation can be clas- 
sified into three types: (1) Individual neurons are oscil- 
lators with varied frequencies, and collective oscillation 
emerges through phase interaction between individual os- 
cillators; (2) Neurons are simple relaxation elements, and 
oscillation emerges as a result of interaction between ex- 
citatory and inhibitory groups of neurons; (3) Neurons 
are instantaneous integrators of incoming signals, and 
oscillation emerges due to signal transmission delay be- 
tween elements. The first two scenarios have been stud- 
ied in detail for many decades, and are already estab- 
lished as fundamental nonlinear dynamical phenomena 
[11) [3, Hj 0, B @1 ■ Light was shed on the third scenario 
rather recently and much is not known mathematically 
in spite of the ubiquity of signal transmission delays, not 
only in networks of biological neurons but also in net- 
works of artificial electrocircuits 0, H, Q • 

In biological neural networks, the transmission delay 
is a sum of axonal, synaptic and dendritic delays. It has 
been reported that the delay can be comparable to or 
longer than somatic membrane time scale [l^, [HI fl^ . 
We consider here a model network in which the trans- 
mission delays are distributed in a large interval of time. 
We model the neuron for simplicity as a discrete-time 
threshold element that updates its state at each time step 
according to the summed input signals, each of which has 
arrived with a fixed transmission delay. 

In the studies of neural networks, synchronous update 
from the preceding states of other elements has often 
been discussed [H Hi, M, M, \l^M,\M,M,^- This 
synchronous update rule can give rise to a period-two 
macroscopic oscillation as the strength of the inhibitory 
connections is increased [Tsl . [II [S [Hi- In the present 
study, we consider the more general case in which inter- 
neuronal transmission delays are distributed widely in 
time. We derive a recurrence equation of the macroscopic 
order parameter representing the mean activity from the 



dynamics of the individual threshold elements. 

The network is found to exhibit a collective oscillation 
with period close to the mean transmission delay. The 
collective dynamics can not only be a simple alternation 
of the consecutive firing and consecutive resting, but also 
a nontrivially sequenced series of firing and resting, re- 
peating in a given period of time. For the system whose 
intcr-ncuronal transmission delays are distributed uni- 
formly in a given range of time, we obtain multiple stable 
periodic orbits. Due to the multi-stability, one can store 
a nontrivial firing sequence in the network. For the case 
of non-uniformly distributed delays, we solve the recur- 
rence equation and find that the network can exhibit not 
only periodic but also quasiperiodic or chaotic dynamics. 

The present paper is organized as follows: In section 
II a recurrence equation of the mean activity is derived. 
In section III the stationary solution for the macroscopic 
dynamical equation and its linear stability is examined 
analytically. In section IV the macroscopic state equa- 
tion is solved numerically for a variety of distributions of 
the delays. In section V the numerical simulation of the 
original microscopic equations is carried out and com- 
pared with the solution of the macroscopic state equa- 
tion. In section VI we discuss the significance of delayed 
networks. 



II. DERIVATION OF A MACROSCOPIC STATE 
EQUATION 

In this section, we derive the dynamical equation of 
the macroscopic activity from the microscopic dynam- 
ics of the individual elements, each of which is a simple 
threshold neuron that evokes an all-or-nothing response 
to an input Vi(t) at discrete times t = 0, 1, 2, 3, • • ■ as, 

Xi{t) = sgn{vi{t)) , (1) 

where sgn(t;) is the sign function that takes values -1-1, 
0, and —1 respectively for v > 0, v — 0, and v < 0. 
Here, Vi(t) is a "total input" representing a linear sum of 
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incoming signals from other neurons, 

n 

where j is the "synaptic weight" that represents exci- 
tatory (positive) or inhibitory (negative) connection from 
the jth neuron to the ith neuron. dij{= 1,2,3- ••) is 
the signal transmission delay specified for each transmis- 
sion line (FigllJ. Si will be called the "external stim- 
ulus" in this paper. In an asymmetrical representation 
in which neuronal firing and resting are represented as 
{xj -I- l)/2 = 1 or 0, the external stimulus should be read 
as Si - Yl%i '^i-r 




time 

FIG. 1: Every neuron receives signals that arrive after de- 
lays in the individual transmission lines. The response signal 
is sent back to the other neurons with different delays and 
different connection weights. 

A macroscopic order parameter representing the mean 
activity is defined as the average neuronal state at each 
time t: 

1 " 

X{t)^-Y,Mt). (3) 

i=l 

In the following, we derive a dynamical equation of the 
macroscopic state, in parallel with Amari's derivation for 
the synchronous update rule p^ . which corresponds to 
a particular case of the unit-time delay in the present 
model, {dij = 1}. 

The macroscopic state defined above is identical to the 
difference of ratios of positive and negative total inputs Vi 
to individual neurons (i = 1, 2, • • • , N). Using the distri- 
bution pt{v) of total inputs v at time t, the macroscopic 
state can be represented as 

/■oo />0 

X{t) = / pt{v)dv~ / pt{v)dv. (4) 

JQ J -oo 

The central limit theorem holds for the summed inputs, 
Wi^jXj {t — dij)}i, in the limit of a large number of 
neurons, if the individual signals {wijXj{t — dij)}ij are 
independently sampled from a given distribution of a fi- 
nite variance. In the present deterministic model, the sta- 
tistical independence holds if individual neuronal states 



{xj } are chosen independently from synaptic connections 
{wij}. As the microscopic states are updated through 
synaptic connections, however, {xj} are not independent 
from {wij} any more as time goes by. Even with this con- 
dition, there would be room for statistical independence, 
if the connection weights {wi_j}i_j themselves are chosen 
independently of each other. The question of statistical 
independence was initially raised by Rozonoer [22]. It 
was proven by Amari that the statistical independence 
holds under some special conditions [1^. In the present 
paper, we further introduce the distribution of transmis- 
sion delays {dij}. This raises another problem of the 
correlation between {xj} and {dij}. In this paper, we 
do not go into this open problem, but rather, use the 
assumption of their statistical independence to construct 
a macroscopic state equation for our proposed delayed 
networks. We will examine the validity of the assump- 
tion by comparing the solutions of the macroscopic state 
equation with the simulation of microscopic equations. 

In addition, if the "external stimuli" {si}i are normally 
distributed, then the distribution pt{v) of total inputs 
V is Gaussian, characterized solely by the mean /it and 
variance at each time t as 

By inserting this into Eq. ^ , the macroscopic state equa- 
tion is obtained as 

Xit) = ^ (^) ' (5) 
where F{x) is the error function: 

F{x) = e-'^dx. (6) 

Under the above-mentioned assumption of statistical 
independence, the mean fj,t and variance of input sig- 
nals V are obtained as 

fit — nwd{t) + s, (7) 
a? = nw^l - d{t f) + nal + al (8) 

where w, s, cr^ and are the means and variances of 
{wi.j} and {sj}. d(t) denotes the mean past activity: 

^ n n 

"W^^^EE^^^^-'^'j)- (9) 

i=i j=i 

We consider the case that delays {dij} are randomly 
distributed from d 1 to m over transmission lines with 
the ratios Pd > 0, {pi+ P2 + ■ ■ ■ + Pm = Assuming the 
statistical independence between individual states {xj} 
and delays {di.j}, the mean past activity is given by the 
weighted average of the past macroscopic states: 

ni 

d{t) = Y,PdX{t-d). (10) 

d=l 
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The evolution equation of the macroscopic state is 
given by inserting the mean //^ and variance CTj into 
Eq.®. If the model parameters satisfy nw^ <C nal^+a^, 
the evolution equation simplifies to 

Xit)=F(wJ2pdX{t~d) + Sy (11) 

where W and S are dimensionless parameters respec- 
tively representing the average synaptic weight and the 
external stimulus, 

W = nw/y/nal+al (12) 
S = s/^nal + al (13) 

We will analyze the recurrence equation (jlip in the fol- 
lowing sections. It should be noted that nu>^ <^ na'^ + a'l 
is not an essential condition for a macroscopic equation 
^ to hold but is merely introduced to make the analysis 
simpler. 

III. LINEAR STABILITY ANALYSIS OF 
MACROSCOPIC STATIONARY STATES 

Given a macroscopic stationary state X{t) = Xq that 
satisfies 



The phase space of parameters W-S is divided into three 
regions according to qualitative differences in the macro- 
scopic state stability (Fig[21); [monostable region]: The 
system has only one stable stationary state (|/3| < 1); 
[bistable region]: The system has one unstable station- 
ary state with /? > 1 and two stable stationary states 
(1/3| < 1); [oscillatory region]: The system has one unsta- 
ble stationary state with f3 < —1, and one stable oscilla- 
tory orbit of period two. 
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FIG. 2: Phase diagram for the synchronous update rule. W 
and S are dimensionless parameters respectively representing 
the average synaptic weight and external stimulus. See the 
text for details. 



X„ = F {WXo + S) , (14) 

we wish to analyze its stability. For this purpose, the 
recurrence equation ()11|) is linearized with respect to the 
deviation from the stationary state, SX{t) = X{t) — Xq, 
as 

m 

6X{t)=f3j2pdSX{t-d), (15) 

where 

X=Xo 

The stationary state is locally stable if all roots of the 
characteristic equation, 

a'"-/3(pi«™-i+P2a'""' + --- + p™-i« + Pm) = 0, (17) 

are smaller than 1 in magnitude. 

A. Synchronous update rule, or unit-time delay 

We start with the synchronous update rule originally 
studied by Amari [1^, which corresponds to the case 
of unit time delay {dij = 1} in the present framework. 
The linearized equation for this case is simply given by 
X{t) — (3X{t — 1). The stationary solution is stable if 
< 1, marginal if |/3| = 1 and unstable if \(3\ > 1. 



B. Uniformly distributed delays 

Next, we consider the case that delays are uniformly 
distributed in a given interval of time, pd — l/m for 
d = 1, 2, • • • ,771. The characteristic equation for this case 
is 

a™- A(a'"-i + a™-2 + ...-^a + l) =0. (18) 

In the following, we prove that for this particular case 
the stationary solution is locally stable if —7n </?< 1, 
and unstable otherwise. 

For /3 — —771, all roots of the characteristic polyno- 
mial equation (jlSp align on a unit circle in a complex 
plane; a = ^^^^k|(ra+\■) ^ 1,2,- •• ,m). As /3 passes 
—m from above, all the roots simultaneously leave the 
unit circle (Fig[3Ua)). For (3—1, the characteristic equa- 
tion has one root a = 1, which has the largest length. 
As P passes 1 from below, this root exceeds 1 along the 
real axis (Figl3]^b)). In Appendix A, we prove that the 
characteristic equation p8|) possesses no other roots of 
unit length. This means that the linearized equation can 
be destabilized only at these two critical points /3 = — tti 
and (3 — 1. 

As in the case of the synchronous update rule or unit- 
time delay {dij = 1}, the phase space can be divided 
into three characteristic regions according to qualita- 
tive differences in the linear stability. The distributed 
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FIG. 3: Roots of the characteristic polynomial equation in a 
complex plane for the case of uniformly distributed delays in 
an interval of m = 9. (a) All the roots leave the unit circle 
as j3 passes — m from above. Filled circles, crosses and filled 
squares represent roots of characteristic equations of stable, 
marginal and unstable cases, respectively, (b) One root ex- 
ceeds 1 along the real axis, as (3 passes 1 from below. 



inside the unit circle in the complex plane. At /3 = 1, one 
root arrives at a = 1, above which a stationary solution 
becomes unstable. 

The critical value of f3 on the negative side of W de- 
pends on the distribution of the delays {pd}d- As de- 
picted by FigOJa), all roots do not necessarily leave 
the unit circle simultaneously, unlike the uniformly dis- 
tributed case. 
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FIG. 4: The linear stability boundary between monostable 
region and oscillatory region shifts to the lower value of W as 
the distribution interval of the delays m is increased. 



transmission delay doesn't shift the boundary between 
monostable and bistable regions. The boundary between 
monostable and oscillatory regions is shifted to the lower 
direction in ly as m is increased (Fig 2]). Note that the 
present categorization is solely based on the (local) linear 
stability and there is room for other dynamical orbits to 
coexist globally. The coexistence of multiple orbits will 
be discussed in the next section. 



C. Non-uniformly distributed delays 

Finally, we consider the general case that delays are 
distributed non-uniformly in a given interval of time m, 
satisfying pi+p2 + ' ■ ■+Pm = 1- Even in this general case, 
the boundary between monostable and bistable regions 
is the same as the case of uniformly-distributed delays as 
well as the unit-time delay, as proven below: At /3 = 0, 
all roots of the characteristic equation (fT7)) degenerate at 
a = 0. While for — 1 < /? < 1, the characteristic equation 
does not have roots of |a| = 1, implying that all roots are 



FIG. 5: Roots of the characteristic polynomial equation, of an 
example of nonuniformly distributed delays, pj — j/45, j = 
1,2, ■•• ,9. (a) Two complex conjugate roots cross the unit 
circle while others remain inside, as /3 passes some negative 
critical value from above, (b) One root exceeds 1 along the 
real axis, as /? passes 1 from below. 



IV. NUMERICAL ANALYSIS OF THE 
MACROSCOPIC STATE EQUATION 

In this section, we solve the recurrence equation (fTTj) 
to observe the dynamics of the macroscopic order pa- 
rameter X(t). The dynamical state obtained for the syn- 
chronous update model or the unit-time delay is a simple 
period-two oscillation. Special attention is paid here to 
nontrivial temporal activity patterns of the network with 
distributed delays. 

A. Uniformly distributed delays 

First, we numerically solve the macroscopic state equa- 
tion of the network with uniformly distributed delays, 
Pd = l/m for d = 1, 2, ■ • ■ , m: 

X{t)^F(w ^^-'^^''^K sy (19) 

Oscillation occurs in the parameter region where in- 
hibition dominates as we see in FigsH] and 21 Figure 
[5] displays the bifurcation diagram of the macroscopic 
recurrence equation in the case of m — 6. For a 
given mean connection that is negative, (a) W = —10 
or (b) W = —20, we vary the external stimulus S and 
observe the temporal activity pattern generated by the 
recurrence equation: For each set of parameters of W and 
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we choose 100 random initial conditions, and iterate 
the recurrence equation for t — 10000 and plot the last 
several values of X{t). 

If the external stimulus S is sufRciently small, the sys- 
tem exhibits a stationary sequence of negative X close 
to —1. As S is increased, a positive X close to +1 ap- 
pears among Xs close to —1, once every m + 1 iterations. 
Note that the oscillation is observed already in the pa- 
rameter region in which the stationary state is linearly 
stable. In other words, oscillatory orbits and a stationary 
state coexist in the same system, as mentioned in section 
III. As S is increased further, a positive X appears twice 
every m + 1 iterations. The temporal order of positive 
and negative Xs can be chosen arbitrarily, but it is fixed 
once the iteration starts. The proportion of positive Xs 
increases with 5*, and eventually the system exhibits a 
stationary sequence of positive X. Note that different 
proportions of positive and negative Xs may coexist for 
the same value of S. 

We solved the recurrence equation from initial condi- 
tions with Xs uniformly distributed. Within our range 
of numerical investigation, we obtained all possible ar- 
rangements of positive and negative Xs of cycle length 
m+1. In the case of m = 6, the period of oscillation is 
always m+1 = 7, the prime number. In the general case 
of m -|- 1 that is not a prime number, the period could be 
shorter with the repetition of short sequenced activity. 
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FIG. 6: The bifurcation diagrams of the macroscopic recur- 
rence equation (|19|) displaying the order parameter X{t) for 
each value of external stimulus S. (a): W = — 10, m = 6; (b) 
W = —20, m = 6. For each value of S, we take 100 random 
initial conditions and plot the last seven X{t). The linear 
stability boundaries for (a) and (b) are Sc = ±6.2 and ±18.2, 
respectively depicted as dashed lines. 



B. Non-uniformly distributed delays 

Next, we examine the case that the delays pd are un- 
evenly distributed: 



x{t) = F ( w ^''=^,^ ^' \ — ^ ± s 



(20) 



where is drawn from a Gaussian distribution with 
mean zero and variance e^. 
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FIG. 7: Top: The bifurcation diagram of the macroscopic 
equation (|2Q[) with bifurcation parameter S. Bottom: the 
Lyapunov exponents, (a) and (b): m — 6,W — —10, e = 0.05; 
(c) and (d): m = 6,W = —20, e = 0.1. For each value of S, 
we take 100 random initial conditions and plot the last seven 
X{t) for each. 



In addition to periodic oscillation, the system may ex- 
hibit quasi-periodic dynamics (FigslZl^a) and (b)). As the 
deviation e is increased further, the system may exhibit 
chaos (FigslZl^c) and (d)) characterized by the positive 
value of the Lyapunov exponent: 



1 |^X(t)| 
A = lim lim - log , , 
i^oo5x(o)^oi |'5X(0) 



(21) 



where SX{t) = {SX{t + l),SX{t + 2),- ■ ■ ,SX{t + m)) is an 
m-dimensional perturbation vector added to an original 
orbit. 



Figure EKb) shows the solutions for W = —20. It is 
notable that there appears a new intermediate state of 
X close to in addition to states of X close to +1 and 
— 1. In comparison with the case W = —10, the re- 
gion in which different proportions of positive and nega- 
tive Xs coexist becomes relatively narrower. In the limit 
W ^ —oo, the proportion of positive and negative Xs 
are uniquely determined by the parameter S' = mS/\W\, 
as is proven in Appendix B. 



V. NUMERICAL SIMULATION OF 
MICROSCOPIC DYNAMICS 

In this section, we carry out a numerical simulation 
of the original microscopic equations ([1]). The dynamics 
of macroscopic order parameter representing the mean 
activity of microscopic elements are compared with solu- 
tions of the macroscopic state equation pTjl . 

The size of the network we examined is n = 1000. 
Synaptic connections {wi^} are drawn independently 
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FIG. 8: Left: The mean activity X obtained by the simulation 
of the microscopic equation: (a) and (b): w — —0.08, aJJ, — 
0.09, s = 0, and al = 0; (c): w = -0.12, al = 0.09, s = 0, 
and o"s = 0. Right: Solutions of the macroscopic equation 
exhibiting similar temporal patterns: (d) and (e): W = —8.4 
and 5 = 0; (f): W = -12.6 and S = 0. The dashed squares 
in (c) and (f) are magnified in Figs[9ja) and (e). 



from a Gaussian distribution of mean w and variance cr^ . 
In this simulation, external stimuli {si} are taken as 0. 
The transmission delays {di.j} are chosen randomly from 
{1, 2, • • • , 6}. This corresponds to the case of uniformly 
distributed delays, pi — P2 = ■ ■ ■ = P& = ^ / ^ 

The resulting macroscopic order parameter X{t) is de- 
picted in FiglSl (a) and (b) display the mean activity in 
the case of w = -0.08, crj = 0.09, s = 0, and a1 = 0. 
The macroscopic state equation (jlip with the parameters 
W = -8.4 and 5 = estimated with Eqs.n!]) and (fTa]) 
exhibits temporal mean activity patterns similar to the 
simulation results (FigslSJd) and (e)). It is interesting to 
observe that the temporal activity pattern is gradually 
modified as time goes by (from (a) to (b)). This would 
be due to the finite size effect. The oscillation is sta- 
bilized in the parameter region of smaller ui, even in a 
system of the same size. The simulation result with the 
parameter w = —0.12 and a solution of the macroscopic 
equation with the corresponding parameter W — —12.6 
are displayed in Figsl5]Jc) and (f). 

As is predicted in the preceding section, the macro- 
scopic order parameter exhibits a wide variety of tem- 
poral dynamics depending on the initial condition. 
FigsiHa)-(d) depict examples of temporal sequences re- 
alized in the same network. Note that these sequenced 
firing patterns are stable, once they arise from a given ini- 




FIG. 9: Various temporal patterns of mean activity X exhib- 
ited by identical networks: Left: (a) to (d): Numerical simu- 
lation of the microscopic equation, with the same parameters 
of FiglSjc). Right: (e) to (h): Solutions of the macroscopic 
equation, exhibiting temporal activity patterns similar to the 
simulation results. 



tial condition. FigslS]^e)-(h) display the similar temporal 
sequences exhibited by the macroscopic state equation 

(ED. 



VI. DISCUSSION 

In the present study, we have demonstrated collective 
dynamics exhibited by the neural network whose intra- 
transmission delays are widely distributed in time. In 
the case that delays are distributed uniformly in time, 
the system is found to exhibit collective oscillation with 
a period close to the average transmission delays. The 
network possesses multiple stable orbits of nontrivially 
sequenced series of firing and resting, reverberating in a 
certain period of time. The multistability has also been 
reported for continuous time dynamical systems that are 
accompanied with delay [l^. Ha, [H, [131 ■ In addition to 
the periodic motion, we also found that the dynamics 
can be made quasiperiodic or chaotic by changing the 
distribution of delays. 

The analysis of delayed networks has mainly been con- 
fined to their stability [H, [H, HO, [M HI ■ Recently, the 
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dynamical aspects of the delayed systems are drawin g at - 
tention; some possible relevance in biology [ssl [3^. Tssj: 
controlling systems by utilizing delays in engineering 
(sgI [stI . [ssI . [39| . Transmission delays provide networks 
with potential applications. In order to control them, it 
is necessary to comprehend full aspects of their dynami- 
cal characteristics. 



Appendix A 

We prove that the characteristic polynomial equation 

m 

_P_ / m-1 ^ + . . . + a + l) = 

TO 

has no roots of unit length, other than the roots a = 
g27r»fc/(m+i) (fc = 1^ 2, • ■ • , to) at /3 = -TO and a = 1 at 

/3 = 1. ^ 

It is readily seen that a = I can be a solution of the 
characteristic equation, if /3 = 1. If a 7^ 1, the character- 
istic equation can be transformed into 

a™+i -1 = + (22) 
m 

This equation means that 1, a™ and 0;™+^ are aligned 
on a line in the complex plane. In addition, if a (a™, 
Q^+i^ lies on the unit circle, at lease two of 1, a™ and 
^m+i jjjy^gi; Idcntical. This is satisfied with Eq. (|22|) 
only if = 1 (a ^2^ik/{,n+i) ^ = 1, 2, ■ • • , to) and 

(3 = — m. 



where S' = mS/\W\. We obtain exact solutions of this 
recurrence equation. For simplicity's sake, we consider 
here the case of noninteger 5", with which X{t) takes the 
value of either -fl or —1. 

The recurrence equation generates a stationary se- 
quence of if S" > to, and a sequence of —1 if S" < —to. 
A sequence composed of both +1 and —1 is generated if 
\S'\ < TO. We prove here that the recurrence equation 
generates a sequence composed of -l-ls and —Is with a pe- 
riod of TO-l-1, with the number of -l-ls being [(to-|-S") /2] : 
This sequence of length to -I- 1 satisfies the relation 

m 

^X{t - d) = -m - 1 + 2\{m + S')/2l 
This can be rewritten as 

m 

X{t)^ - ^X{t-d)-m-l + 2\{m + S') /2] . (24) 

d=l 

For X = ±1, this equation is identical to Eq. (|^51) . This 
means that an arbitrary sequence of a period of to -I- 1, 
with the number of -l-ls being [(to -|- S")/2] is a solution 
of the recurrence equation 
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